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Abstract. In this paper, we propose an efficient implementation of combining 
Dynamical Mean field theory (DMFT) with electronic structure calculation based on 
the local density approximation (LDA). The pseudo-potential-plane- wave method is 
used in the LDA part, which makes it possible to be applied to large systems. The full 
loop self consistency of the charge density has been reached in our implementation 
which allows us to compute the total energy related properties. The procedure 
of LDA+DMFT is introduced in detail with a complete flow chart. We have also 
applied our code to study the electronic structure of several typical strong correlated 
materials, including Cerium, Americium and NiO. Our results fit quite well with both 
the experimental data and previous studies. 
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1. Introduction 

The first principle calculation based on the density functional theory (DFT) with 
the local density approximation (LDA) and its generalization generalized gradient 
approximation (GGA) is very successful in predicting ground-state properties and band 
structures of a wide range of real materials. However it is known for a long time 
that it is not sufficient to calculate the electronic structure of those strongly correlated 
materials (i.e. transition metal oxides, actinide and lanthanide-based materials, high 
Tc superconductors) by these methods alone even on the qualitative level. In order 
to overcome these shortcomings of the traditional DFT-LDA scheme, remedies such as 
LDA+U have been proposed, which can describe some of the strong correlated materials 
with long range order in their ground states. While even LDA+U can not be applied to 
many of the strongly correlated materials, for example, the paramagnetic Mott insulator 
phase and the materials containing unfilled f-shell with strong multiplet effects. 

Therefore it is important to have a first principle method which can be applied to 
the strongly correlated materials with the ability to capture the dynamical properties, 
finite temperature properties and multiplet effects. In the last two decades, the 
dynamical mean field theory (DMFT) has been developed very fast to be the standard 
tool to deal with the on-site correlation effect in the limit of large dimension [1]. 
After being successfully applied to many model systems, DMFT has been considered 
as a powerful tool to capture the on-site correlation effects based on the Hubbard 
like Hamiltonians containing both the local interaction terms and the single particle 
Hamiltonians extracted from LDA. Therefore the LDA-I-DMFT method which combines 
the DFT-based band structure techniques with DMFT has been proposed and developed 
quickly in the last decade. By DMFT the local correlation effect can be well described by 
the self energy, which has frequency dependence and in general takes the matrix form 
within the subspace spanned by the correlation orbitals. Using the Green's function 
containing self energy, many of the physical properties of strongly correlated materials 
can be calculated, i.e. the electronic spectral function, the total energy, the optical 
conductivity and the local spin susceptibility. 

Most of the LDA-I-DMFT calculation till now have been performed by the partially 
self-consistent scheme, where the local self-energy is obtained by the DMFT calculation 
with a fixed LDA Hamiltonian generated from the fixed LDA charge density. Therefore, 
in this simplified scheme one neglects the inference of the strong on-site Coulomb 
interaction on the charge density. The above mentioned "one-shot" DMFT calculation 
works quite well for the electronic structure. While in order to obtain reliable total 
energy related properties, i.e. the equilibrium volume, the elastic constants and the 
phonon frequencies, the LDA-I-DMFT scheme with full charge density self consistency 
is needed. Up to date, there have been several fully self-consistent LDA-I-DMFT schemes 
[2, 3] as well as actual implementations [4, 5, 6, 7]. As mentioned in [7], there are three 
major issues that have to be addressed in DFT-I-DMFT implementations: i) quality of 
the basis set, ii) quality of the impurity solvers, iii) choice of correlated orbitals onto 
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which the full Green's function is projected. In the implementations of LDA+DMFT 
mentioned above with full charge self-consistency , the basis set of the linear muffin-tin 
orbitals (LMTO) are used, and the local correlated orbitals are naturally chosen based 
on the muffin-tin orbitals . While on the other hand, there are very few implementations 
of LDA+DMFT with full charge self-consistency are based on the plane wave methods. 

In the present paper, we implement a full self-consistent LDA+DMFT scheme 
using pseudo-potential-plane- wave as the basis set, and we use either atomic orbitals or 
Wannier orbitals depending on different physics systems. The Hubbard-I approximation 
is used as the impurity solver in the present paper, but it can be replace by more precise 
solver like continuous time quantum Monte Carlo (CTQMC) for metallic systems. The 
present paper is organized as follows. In section 2, we describe the way to choose and 
construct the correlated orbitals onto which the local interactions exert. In section 3, 
we show the full- loop flow and LDA+DMFT (Hubbard-I) formalism in detail. We apply 
this LDA+DMFT approach to several correlated materials such as Ce, Am, and NiO in 
section 4. Finally we conclude our work in section 5. 

2. Projection onto localized orbitals 

In the LDA+DMFT calculation, all the energy bands are divided into two groups. 
And only the on-site interactions among those local orbitals are treated more precisely 
by DMFT. Therefore the choice of localized orbitals does affect the results obtained 
by LDA+DMFT method and becomes one of the important issues in LDA+DMFT. 
A natural choice for the local orbitals is a set of atomic like wave functions with 
corresponding d or f characters. Typical atomic like local orbitals is the linear muffin- 
tin orbitals (LMT0s)[8] adopted in early LDA+DMFT implementations [9, 10, 11]. 
However a more physical choice should be wannier functions, since the shape of the local 
orbitals will be altered in crystals especially when there is strong hybridization between 
localized orbital and delocalized p-type or s-type orbitals. The wannier functions are 
not uniquely determined by the Bloch wave functions, so different choice can be made, 
for example, the Nth-order muffin-tin orbitals[12], the projected wannier functions[13] 
and the Maximally localized Wannier functions (MLWFs)[14, 15]. Comparisons of these 
choices have been made in previous literature. In this paper, two kinds of local orbitals, 
atomic like functions and the projected wannier functions are adopted according to the 
different systems. 

This implementation of LDA+DMFT method reported in this paper is developed 
on an existing DFT package ESTATE (Beijing Simulation Tool of Atomic TEchnique), 
which is based on the pseudo-potentials and plane waves. Unlike the LMTO methods, 
the local orbitals do not enter the basis set of plane waves and should be constructed 
in a suitable way. Projected wannier functions or MLWFs as local orbitals have been 
used by previous reported LDA+DMFT implementations based on pseudo-potentials 
(or projector augmented wave (PAW)) method) and plane wave scheme[16, 17, 18]. In 
this report, two kinds of orbitals are used according to the need, one is directly derived 
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from the atomic wave function of an isolated atom and the other is projected wannier 
functions constructed from these atomic wave functions. In order to be self-contained 
the construction procedure is presented below in detail. Since all our calculations are 
based on plane wave set, the local orbitals will be projected onto these plane waves. 

First we consider the atomic like local orbitals. The localized nature of the 
correlated bands with d or f characters in crystals ensures that these local orbitals 
are very similar to the corresponding atomic wave functions of an isolated atom. The 
atomic wave functions can be picked as the local orbitals directly if the hybridization 
in crystal could be neglected. All atomic wave functions for the isolated atom could be 
obtained by the solving the all-electron radical Schrodinger equation 

{T + V{r))^^,'^^{r) = E^,^^l^^{r) (1) 

In this manner, all orbitals are well defined by the primer quantum number n and 
angular quantum number /. In realistic systems, the local orbitals usually come from 
the partially filled dor f shell and could be picked according to its character. Often there 
is only one type of local orbitals, so the local subspace can be labeled by /. Of course, 
considering the spherical symmetry of the isolated atom, the spherical harmonics should 
be multiplied the radical wave functions to form a complete set. The local orbitals are 

\M = \<j''yim). (2) 

and in which m = 1 . . . 2/ + 1 denote different angular components. 

In methods based on plane waves, it is not desirable to use the all-electron wave 
function directly since it requires a large amount of plane waves to expand in momentum 
space. To avoid this problem, in PAW method all-electron atomic partial waves can 
be used, while in pseudo-potentials method, pseudo atomic wave functions can be 
chosen. The latter is used in this paper. During the generation of pseudo-potentials, 
the Schrodinger equation obeyed by the pseudo wave functions is as below 

iT + V^"(r))V^„7(r) = E^Jij^ir) (3) 

The details on the generation of pseudo-potentials can be referred to previous 
literature[19, 20]. The spirit of pseudo-potentials is quite straight forward. Beyond 
a core radius r^, the pseudo-potential V^^ play the role of a scattering center just as a 
real atomic potential V^^^ do. Thus, the pseudo eigen energy E^^ should be the same 
as the realistic one E:^/"^, and both the pseudo wave functions iV'nf) the pseudo- 
potential V^^ coincide with the exact wave functions \ipn^^) realistic potential 
yALL q£ isolated atoms beyond this core radius Tc, respectively. 

The quality of pseudo wave functions is the same as the quality of the pseudo-potentials, 
which could be justified by comparing simple DFT calculations with accepted results. 
The pseudo wave functions bear the same atomic features as the exact atomic wave 
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functions, and can be picked as the local orbitals as above (angular quantum number / 
ignored since usually only one kind of local orbital are considered). 

It is convenient to transform the local orbitals from real space to momentum space in 
crystal calculations 

\(t>rn,i.)=Y.^''''''\^r,^,I)■ (6) 
/ 

Since the atomic wave functions at different site / are not orthogonal, an 
orthogonalization procedure is needed for the above local orbitals. 

The DFT calculations give a set of Bloch waves |\E'nk) spanning the total Hilbert 
space in which lies the local subspace spanning by the local orbitals. Thus, the physical 
choice of local orbitals is Wannier functions derived from the Bloch waves. The Wannier 
functions are localized in real space and could be constructed from the Bloch waves via 
a unitary transformation. However, this unitary transformation is not unique since the 
phase factors of Bloch waves are uncertain, which results in that there are many kinds of 
Wannier functions in use, e.g., the projected Wannier functions, NMTOs and MLWFs. 
Also it is not necessary to include all the Bloch waves but the relevant bands lying in 
a certain energy range to construct the Wannier functions. The local orbitals could be 
picked as a subset of these Wannier functions with specified d or f localized character. 

The projecting Wannier functions method is a simple way to construct wannier 
functions, in which trial wave functions are projecting onto physical relevant Bloch 
bands. The atomic like local orbitals above are used as trial functions here. The 
Bloch bands can be selected by band indices {Ni, ...Nj) or by an energy interval {Ei, Ej) 
enclosing them. 

|W^m,k) = IV^nk)(V'nk|0«m,k) 
n=Ni 

= l^nk)(^nk|0m,k)- (7) 

n(Ei<e„k<£^j) 

These orbitals need normalization and orthogonalization to be true Wannier functions. 

|W^n^,k) = ECm,m'(k)~'/Wm,k) (8) 
m' 

in which O is the overlap matrix between different |PFm,k)s 

Om,m'{k) = {WmM\Wm',^) (9) 

Although the construction method can not give the most localized Wannier orbitals, the 
basic features of the localized bands are captured to a very good extend as indicated by 
a few reports[13, 16, 17, 18], and also proved by our calculations in this paper. 
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3. LDA+DMFT(Hubbard-I) Formalism 

We will introduce the detail procedge of LDA+DMFT with full loop charge density self 
consistency in this section. First we plot the general flow chart in figure 3, which contains 
a inner DMFT loop and a outer charge density loop. However, if we use Hubbard-I 
approximation as the impurity solver, the inner DMFT loop is thus neglected. 

In the following subsections, we are going to describe the whole process in detail. 

3.1. First LD A Calculation 

The first step of the full loop LDA+DMFT is an LDA self-consistent calculation, whose 
main purpose is to generate an effective single particle Hamiltonian Hlda and construct 
the correlated orbitals. 

Generally the effective LDA Hamiltonian can be expressed as the following: 

Hlda = ^nkclk^nk (10) 

nk 

In the above equation, n = 1 ~ Nf,and are the joint indices of band and spin; Enk is the 
eigen energy of the Kohn-Sham equation determined by HLDAl'^nk) = -E'nkl'/'nk), where 
Iv'nk) represents a set of orthonormal Bloch functions. 

The LDA calculation also complement the construction of local correlated orbitals. 
In the present implementation the atomic orbitals and the Wannier functions are two 
types of commonly used local basis. In LDA+DMFT, the local interactions have been 
considered on the DMFT level only within the local orbitals and in order to set up 
the DMFT self consistent equation we need to obtain the overlap matrix between local 
orbitals and Bloch wave functions, which takes the form of 5"^^ = (oklv^nk); with Greek 
letter a denote the index of local orbitals. The completeness of the Bloch basis set gives 

EKnP = l (11) 

n 

for any local orbital a and any k-point. 

3.2. DMFT Loop 

The purpose of the DMFT loop is to calculate the self-energy caused by the local 
interactions through the DMFT self consistent loop. As we have mentioned before, the 
DMFT loop is not necessity if we use Hubbard-I approximation as the impurity solver. 
However, here we first introduce the algorithm with more general impurity solvers. 

By adding the local interaction terms to LDA, we get the total Hamiltonian of the 
system which is to be solved by DMFT, 

H LDA+DMFT = HldA + ^{^1 — V^c) ~ (12) 

i 

where i is the site index. For a specified site, we remove the i index and 

Hu = \t. u^pJiflf.h (13) 
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first LDA calculation to obtain the LDA effective Hainiltonian 
^LDA = Enkc'lf.Cnk OR the basis set of Bloch functions Itpnk) 

nk 

and to construct the local orbitals by evaluating the overlaps 
between local single particle states and Bloch states {xklfnk)- 



reconstruct Hlda on the basis set of new Bloch 
functions and by the LDA calculation with charge 
density p(r) and local orbitals \xk) unchanged. 



add one physically reasonable correlated interaction 
to LDA effective Hamiltoniaii, to get the total corre- 
lated Hamiltonian IJld.a+da/ ft = Hlda + Hu- 



the full Green's function projected onto local orbitals is given 
by Qf{iuj) = T.{<^k\^>^k)—^^z^)i^^>'\^k) and is regardec 
as the Weiss mean field in the dynamical mean field theory 



the hybridization function is given by A{iLu) = iuj + 1^-- [^o(*'^')l ^ 



the Anderson impurity model into which the lattice model 
maps is Fi^p = T^^qJ^Jqc. + E Vlf^'J^A* + ''-t') + 

Y^E'^^Pflfa + Hu subjected to the hybridization function. 



solve the local 
self-energy 
check if it is 
nverged. 



map the local seU'-energy back to the Bloch basis 
set, and calculate the k-dependent lattice Green's 
function: = M .^_^J^_f,^.^^ \'pmk) 



evaluate the k-space density matrbc p^"^ 
fromGpH:pp = jE'^rl^'^) 



Fourier transformation from to real space charge density p{r) 




calculate quantities: 
total energy ... 



Figure 1. The most general flow chart for LDA+DMFT scheme. The first LDA 
calculation gives out the band structure on LDA level, and constructs local orbitals, 
then the full Hamiltonian is established and solved by a DMFT loop, after which the 
charge density can be recalculated by Fourier transforming the k-space density matrix 
to real space. The "non-interacting" Hamiltonian is thus regenerated based on the new 
charge density profile and the calculation will be completed when the self consistency 
has been reached for the full loop. 
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is written as a general form of two-body local interactions, in which the creation 
and annihilation operators of correlated orbitals fl and are associated with Bloch 
operators in term of 

/l = E4k(^nk|«k) (14) 



nk 



fa = '^Cnk{ak\<^nk) (15) 
nk 

The third term of (12) is the double- counting term correspond to the correlated 
energy that has already been considered in LDA calculation at a Hartree-Fock mean 
field level, 

VDC = T.VScflfa (16) 

a 

and this term will be discussed in section 3.2.3. The last term of (12) 

iV = E4cnk (17) 
nk 

is the total number of particle operator, while the chemical potential fi controls the 
occupation number of the unit cell. 

3.2.1. Quantum Impurity Hamiltonian In DMFT, the correlation problem on the 
lattice can be mapped to a quantum impurity model, which contains the same on-site 
interaction and reads, 

Himp = E ^.Aa^c. + E y^clJ^ + h.C.) + E ET'flfa + Hu (18) 
qa qa a 

The inference from the rest of the lattice site besides the one considered in the impurity 
model is simulated by a non-interacting "heat bath", which is described by the first 
term in (18). The second term describes the coupling between the impurity site and the 
heat bath and the rest two terms describe the local interactions. 



3.2.2. Hybridization function and Weiss field The above quantum impurity model can 
be solved by the impurity solver, i.e. Hubbard-I, and the self energy Il{iuj) is then 
obtained, with which we can construct the lattice Green's function by applying the 
same self energy term as, 

[Gfatuce]-' = iu;-H^- t^iiu) + fJ. (19) 
The hybridization function, which characterizes the dynamics of the "heat bath" 
is defined as A{iu)ai3 = ^ap l^T^ ^-iid can be obtained iteratively by the following 

DMFT self consistent equation. 



A(ici;) = io; — Ei^p — S + 



n -1 



E^«^Mice(i^) 



L k 



(20) 



The above equation is obtained by requiring that the local Green's function on the 
lattice should equal to the Green's function of the quantum impurity problem. The 
equations (19) and (20) form a closed self consistent loop, which determines both the 
self energy and the hybridization function iteratively. 
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3.2.3. Double Counting Term In this section we discuss the exphcit expression of £'^"*J' 
as well as the double-counting term in (12). For 3d system under cubic symmetry, the 
local interaction can be simplified as coulomb interaction U and Hund's rule coupling 
J, which can be written as, 

Hu = UY^ nb,tnb,i + (f/ - 2 J) ^ nb,anb',a - J ^ibafiwa (21) 

b b<b' b<b' 

a a (7 

And the double counting energy in this case has been studied by Held et al [21]. 
and can be approximately chosen as 



Edc — 2 ^^lda+dmft(''^lda+dmft 1) (22) 

where 

^ U + 2{M-1){U-2J)-{M-1)J 

2M - 1 ^ ' 

M = 21 + 1 (24) 

and riLOA+DMFT is the total number of electrons on correlated orbitals for a specific atom, 
hence in (10) becomes 



vdc=y: 



dUa 



— U ("^LDA+DMFT r))''^a (25) 

The final expression of is 

ET" = Y.(ak\HLDA\ak) - V^c (26) 

k 

where 

^DC ~ ^('^LDA+DMFT (27) 

The way to remove the form of double-counting energy is not unique, and in fact 
this process needs intuition of physics. The double-counting discussed in earlier sections 
is a very usual way which is also used in YDA+U method. Besides, the following two 
ways are also commonly used 

• to remove self-energy at zero- frequency S(ia;) — S(ia;) — S(0) 

• to remove self-energy at infinity-frequency S(icj) -+ S(ia;) — S(oo) 

Furthermore, the double-counting energy is also regarded as "impurity solver 
dependent". For example, it is reasonable to use an integer to replace ?t.lda+dmft in 
(22) for the Hubbard-I solver [G], and in this paper we also follow their suggestion. 
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3.2.4- Hubbard-I solver The essential point of Hubbard-I solver is to neglect the effect 
of the heat bath and take the atomic self energy as the zeroth order approximation for 
the quantum impurity problem, which can be written as 

In the above equation (5*^*°™ is the Green's function for a single atom, which can be 
expressed as 

r<atom(- ,N (F")rr'(-F" ^)rT(^r + ^rp /„qx 

<^aa' (l^J - — p , p l^yj 

— -C/r' + 

where |r) and |r') are the atomic eigenstates obtained by exact diagonalization of a 
single atom problem, Xr^ = e~^^^' / {J2r ^~'^^^) represents the occupation probability of 
the local configuration |r), and Fpp, = (r|/Q,|r'). 
Hubbard-I approximation is 

S ^ (30) 

In reference [22] , Savrasov et al propose that the above self energy can be written by 
the summation of a set of poles, which greatly simplifies the DMFT calculation, because 
it is not necessary to handle the full frequency dependence of the Green's function. In 
Hubbard-I approximation, it can be proved that both the atomic Green's function and 
self energy are diagonal and can be written in the form of pole expansion if the following 
two necessary conditions are satisfied. 1) The single particle basis used here should be 
the one which diagonalize the local density matrix; 2) The atomic Hamiltonian only 
contains the two-body interaction terms. Obviously by chosen the proper single particle 
basis the atomic Hamiltonian (18) in general will satisfy the above two condition and 
thus can be always written in terms of pole expansion. Unlike the reference [22], where 
they only use a few poles to capture the main features of the self energy, in the present 
implementation, we keep all the poles in the self energy, which makes it more accurate. 
From (29), it is very clear that the atomic green's function has already been expressed in 
terms of poles and the pole expansion of the corresponding self energy can be obtained 
using (28), which is introduced in details in the Appendix. Therefore in general the 
above atomic self energy can be written as 

np 
i=l 

where i labels the number of poles, np is the total number of poles in the self energy 
and Vi is a vector defined in the orbital space describing the distribution of the ith "pole 
states" among the local orbitals. 



3.3. Correction of the Density Matrices and Pole Expansion of the Self-energy 

Once we obtain the converged local self-energy, we are at the point to correct the k- 
dependent lattice Green's functions as well as the density matrices. In general we need to 
do the summation over the Matsubara frequencies, which is time consuming for realistic 
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materials contains lots of bands. But it can be greatly simplified when the self energy 
can be expressed as the summation of poles, which can be written as 



(V^nik|PLDA+DMFT|v^n2k) 

1 1 
= ^T.(Vn,k\-. — -7-— IV'n^k) (32) 

As pointed out firstly in reference [22], when the self energy can be written in terms 
of poles, the full green's function can be expressed as the "physical" part of an enlarged 
"pseudo Hamiltonian" , defined as 

+ -j*) (33) 

where HldaO^) = -E'nkc]^kC„k is the physical Hilbert space, Vi is defined in (31) and 
P = Diag{pi,p2, ....,Pnp)- Therefore the physical green's function can be expressed in 
terms of eigenstate of the pseudo Hamiltonian simply as 

1 



iun - Hlda{^) + Vdc - S(za;„) 



V^nak) 



= (<^nik|^/'fk): ^ X (^fkl<^«2k) (34) 

where iV'fk) a-nd ^i^O^) the eigenstate and eigenvalue of pseudo Hamiltonian 
respectively. 

Then the sum of frequencies in (32) can be performed directly 

(V'nik I PLDA+DMFT \ </2n2k) 

= ^im(V'niklV^fk): ^^5J77T(V^fkl</'n2k) 

P /k - 

= 5Z(¥'"ik|^fk)(V^fkl¥'n2k)^F[^r(k) -/i] (35) 
ik 

in which yU, is exactly the chemical potential in (12), determined by the occupation 
number of electrons Ntot in the unit cell 

nk 

= El(^nk|OlV[i?r(k)-/i] (36) 

ink 

3.4- Evaluation of Real Space Charge Density and Complete of the Full Loop 

After DMFT obtains the corrected density matrix p, the real space charge density will be 
generated based on it by Fourier transformation, and then, the new LDA Hamiltonian 
as well as the new overlap matrices between Bloch states and local orbitals will be 
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recalculated, which closes the full iteration loop for the charge density self consistency. 
The occupation number of the Bloch basis can also be obtained by, 



n 



LDA+DMFT 



E (ttk|</5„jk) (V^mklplv^nak) (v^n.2k|«k) 
k,nin2,a 



(37) 



3. 5. Calculation of Physical Quantities 

3.5.1. Energy Functional Most of the physical quantities we are interested in can be 
calculated once we have reached the charge density self consistency. First we present the 
energy functional of the full loop LDA+DMFT, which can be used to calculate many 
quantities, i.e. the total energy, force and so on. It can be written as, [6] 

E = Elda[p] - {Hks)lda + (Hks) + (Hu) - Edc (38) 

in which -EldaIp] is the expression of the energy within density-functional theory; 
{Hks)lda = ^t^[HldaGlda] is the non-interacting energy at LDA level; {Hks) = 
trlHiDAG] is the non-interacting energy at DMFT level; (Hu) = |tr[EG] is the 
interaction energy caused by local correlation interactions; the last term Edc is the 
double counting energy given before in (22). In above the meaning for "trace" is defined 
as 



1 



*^[^] = ^im(¥'"k|^(iWn)|'^nk) 



(39) 



nk iuJn 



and the integral path surrounds all the energies of occupation states. 

The first and the second term is easily calculated in the LDA framework, and the 
three terms remaining must be evaluated in the DMFT process, explicitly 



{Hks) = ^Enk{(Pnk\p\Vnk) 



(40) 



nk 



and 



(Hu) 



-tr[SG] 

itr[(Go^-G-i)G] 



E (^mk - /^)l(¥'nk|^mk)| V(^^k " P) 



m,nk 



■{Hks) + {Vdc) 



(41) 



while {Hks) is expressed in (40), and {Vdc) is the expected value of Vdc at DMFT 
level 



{VDC)=J:U^' 



n 



LDA+DMFT 



Uiin 



LDA+DMFT 



1 

2'^ 



LDA+DMFT 



(42) 



(43) 



Notice that here emerges the site index i. We should be very careful in the case that the 
number of correlated atoms is larger than one. The final expression of the total energy 
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IS 



E = Elda[p] — {Hks)lda + {Hks) 
1 



+ 



m,nk 



2 



(44) 



However, if we use Hubbard-I as the impurity solver, and remove double- counting 
term by using an integer, the energy functional is given by Haule[7] 

E = Elda[p] — {Hks)lda + {Hks) 



1 

+ 2 



(^^k - /i)|(^nk|^mk)| V(^^k - P)- {Hks) 



m,nk 



(45) 



where Uq refers to the integer used to remove double counting for the ith correlated 
atom. 



3.5.2. Expressions of DOS and PDOS In order to calculate the density of state (DOS) 
or the partial density of state (PDOS) on correlated orbitals, it is necessary to calculate 
the retarded form of the Green's function. By using the virtue of pole expansion method, 
we write directly 



Die) 



-— )Im 

TT 



and the PDOS of orbital a 

, 1 



DJe 



)Im 



mk 

kin e + /^-^mk + i^ 

(«k|V^mk)(^mk|ak) 



{rj 0^ 



TT 



E 



where 



ttklV^mk) = E('^k|V5nk)(V5nk|^mk) 



(r/ 0^ 



(46) 



(47) 



(4^ 



4. Benchmark 



To validate the LDA+DMFT framework based on pseudo-potentials-planewave package 
and pole expansion of self energy, we benchmark our implementation by applying to 
7-cerium, americium and paramagnetic NiO. These three are typical strongly correlated 
systems in which the valence electrons are believed to be on the localized side as reported 
in previous literatures. We use the Hubbard-I method as impurity solver, which is good 
enough to capture the atomic-like features in the Mott insulators, so we expect that 
these systems could be well described by our method. We would emphasize here that in 
our implementation we can replace the Hubbard-I solver by any of the solvers as long 
as the self energy can be written in pole expansion form. 



14 



4.1. Cerium 

The cerium metal attracts lots of research interests for its isostructural volume-collapse 
transition from 7 phase to a phase. The volume change is about 15% during the 
transition, which is possibly driven by the entropy change [23]. Normal LDA calculations 
could not give a correct description of cerium, as show in table 1, the equillibrium 
volume of cerium given by LDA calculations is smaller than the experimental volume of 
a phase. This is due to the fact that in LDA the /-electrons are treated as itinerant and 
the strong correlation effects among them can not be well captured. When /-electron 
is treated as core electrons, we see that the equilibrium volume is larger but close to 
7 phase. These simple LDA results implies that in 7 phase, the /-electrons are more 
localized than in a-phase. Therefore /-electrons in 7 phase is quite close to the localized 
picture like the situation in Mott insulators, which makes it suitable to applying the 
Hubbard-I solver. The a-Ce which is stable in low temperature is a correlated metal, 
which is confirmed by many other LDA-I-DMFT calculations [2-3], and also by our newly 
development LDA-|-Gutzwiller method[24]. 

In our calculation for 7-Ce, we used ultrasoft pseudo-potentional in which 4/ 55, 
5d, 6s orbitals are taken as valence orbitals. The energy cutoff for plane- wave expansion 
is 12.5 Ha for convergence. Calculations are performed with 10 x 10 x 10 k-points. 
For DMFT calculation, the on-site Coulomb interaction U is chosen to be 6.0 eV as 
suggested in previous reports. [25, 23, 26] 

In figure 2 the total energy versus volume for 7-cerium obtained by LDA, LDA+f/ 
and LDA+DMFT are plotted. The equilibrium volume and bulk modulus are also 
shown in table 1. We can find in figure 2 that the equilibrium volume obtained 
by our LDA+DMFT calculation is very close to the experimental data, which shows 
great improvements over LDA . The bulk modulus is also improved but still larger 
than the experimental data, which may imply the contribution for lattice vibration is 
inneglectable as discussed in reference [6]. 

Table 1. Lattice parameter and bulk modulus of 7-cerium according to both 
experiment and calculations 





Volume (A'^) 


Bulk modulus(GPa) 


Experiments 


34.35[27] 


21[28] 


LDA 


22.36 


62.09 


LDA(fcore)[29] 


36.39 


30.2 


LDA+C//PAW[3(J] 


32.0 


34 


LDA+DMFT [6] 


30.10 


48.49 


LDA+DMFT 


33.29 


38.27 



Figure 3 shows the 4/ partial DOS for 7 cerium at the equilibrium volume, the 
position of the two Hubbard bands agree well with the experimental result. 
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Figure 2. Total energy versus volume curves obtained by LDA, LDA+[/[(i], and 
LDA+DMFT. The dashed lines are fitted by Birch EOS. The energy were shifted for 
different curves. The experimental equilibrium volume is from reference [27] 




Figure 3. 4/ pdos for 7 cerium in LDA+DMFT (Hubbard I), XPS and BIS 
experimental data are from reference [31] and [32] 
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4-2. Americium 

Americium, which is widely used in smoke detectors, can be viewed as the "changing 
point" of the actinide series, where the 5/ electrons changes from delocalized to 
localized states[33]. Valence-band ultraviolet photoemission experiment by Naegele 
et al [3 i] shows that the 5/ states are strongly localized. Am undergoes a series of 
structure phase transitions with the increment of pressure, from double hexagonal close 
packed(dhcp, POa/mmc), to face centered cubic(fcc, Fm3m, 6.1 GPa), face centered 
orthorhombic(Fddd, 10.0 GPa), and to primitive orthorhombic structure (Pnma, 16 
PGa)[3 j]. Therefore Am provides an interesting playfield to investigate the relation 
between /-electron localization and structure transition. Here we focus on the fee phase 
and dhcp phase. 

In the calculation, we used ultrasoft pseudo-potentional contains 5/, Qp, 6d, 7s 
orbitals as valence orbitals, the plane-wave expansion kinetic energy cutoff of 12.5 Ha. 
The LDA self-consistent calculations were performed with 10 x 10 x 10 k-points grid for 
fee phase, 9x9x3 k-points grid for dhcp phase. For the volume calculation, we only 
considered density- density interaction, and on-site Coulomb interaction F^^^ = 4.5 eV, 
for the partial DOS calculation, we considered general interaction, and take the atomic 
values = 7.2 eV, F^^) = 4.8 eV, F^^) = 3.6 eV[36]. 

We present the optimized volume and bulk modulus for both dhcp and fee phase 
in table 2 and figure 4. Our LDA+DMFT results are quite close to the experimental 
data and show great improvement over LDA. Also our results are quite consistent with 
the reference [22], in which the full potential LMTO methods are used in the LDA part. 
The equation of state of Am obtained by our LDA-I-DMFT calculation is plotted in 
figure 5, which again shows very good agreement with both the experimental data and 
previous LDA+DMFT calculation based on LMTO method[22]. 

Table 2. Am equilibrium volume and bulk modulus for dhcp and fee phase by different 
methods, together with the experimental results [S-'j, 37]. 





Volume(A^) 


Bulk modulus(GPa) 


experiment 


29.250 


29.9 


GGA(dhcp)[38] 


19.916 


70.0 


LDA(dhcp) 


18.55 


118.99 


LDA(fcc) 


17.73 


169.40 


LDA+DMFT(dhcp) 


27.04 


52.70 


LDA+DMFT(fcc) 


26.99 


50.75 


LDA+DMFT[22] 


27.4 


45 



Figure 6 shows the 5/ partial DOS of fee phase Am obtained both by our 
LDA+DMFT calculation and LDA+DMFT(OCA) by Savrasov[22], together with the 
photoemission data from ref. [31]. Our results is quite consistent with the previous 
calculations and the experimental results. 
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Figure 4. Total energy per atom of Am obtained by LDA and LDA-I-DMFT. The 
LDA-|~DMFT results of Am for both phase are very close, and the experimental 
equilibrium volume is denoted by dot line. 




Figure 5. Calculated P-V relation of dhcp and fee Am. Here, Vo is experimental 
equilibrium volume. The experimental P-V results [3 7] are shown by black dots. 
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4-3. Paramagnetic phase of NiO 

NiO has been heavily studied as a prototype of Mott insulators. It has been studied 
within the frame of LDA+DMFT by many groups[39, 40]. Therefore NiO can be used as 
a good benchmark material for the implementation of LDA+DMFT. In this paper, we 
apply our code to study the paramagnetic phase of NiO. First we plot the partial density 
of states (PDOS) obtained by LDA in figure 7, which clearly shows metallic behavior 
and is not consistent with the experiments. In our LDA+DMFT(Hubl) calculation, we 
use the wannier functions as the local basis and choose on-site interaction U = S.OeV, 
Hund's coupling J = l.OeV. The Mott insulating features of NiO then can be well 
captured by our LDA+DMFT method, as shown by the PDOS in figure 8). 

The energy gap obtained by our calculation is around 4.0eV, which is also quite 
consistent with the experiments. The photo emission and BIS data obtained by 
Sawatzky and Allen[il] are also plotted in figure 9. We can find that our results fit 
the photo emission and BIS data very well. 

5. Conclusions 

The new implementation of LDA+DMFT based on the pseudo-potential plane-wave 
method is introduced in detail in this paper. We choose the Hubbard-I method as the 
impurity solver to solve the quantum impurity problem generated by DMFT, which 
is quite suitable for the Mott insulator materials. The most important advantage of 
Hubbard-I method is that it is simple enough, which makes the full loop charge self 




^ 1 1 — I 1 1 1 1 
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Energy(eV) 



Figure 6. The spectral function of Am evaluated by LDA+DMFT. The experimental 
photoemission results is from reference [■U] are denoted by solid dots. 
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Figure 7. The partial DOS of NiO obtained by LDA. The main contribution to the 
Fermi surface is attributed to eg-likc Wannier orbitals, and the t2g-like orbitals are 
almost fully occupied. 




Figure 8. The spectral function of NiO obtained by LDA+DMFT calculation. As 
shown in the figure, compared with LDA results, the t2g-like orbitals are still fully 
occupied while the eg-like orbitals split into upper and lower Hubbard bands. 




Energy (eV) 



Figure 9. The spectral function of NiO evaluated by LDA+DMFT, with the 
comparison of experimental data[41]. The behavior of the density of state near the 
fermi surface fits well with the experimental data. 



consistent calculation accessible. We also point out in the paper that the difficulty 
of handling frequency dependent Green's function can be completely removed by 
expressing the self energy in terms of pole expansion, which greatly raise the efficiency 
of the method. Finally we benchmark our implementation of LDA+DMFT on several 
important correlation materials including 7-Ce,Am and NiO. Our results for all these 
materials fit very well with both the experimental data and the previous LDA+DMFT 
results. 



Appendix A. Self-energy in pole expansion form 

In this appendix, we will introduce how to evaluate self-energy from Green's function 
in the pole expansion form. 

Green's function can be expressed as 

GM^E-^ (A.1) 

where Pi is the ith pole of the Green's function, Vi is the weight of the pole, and Nq is 
the number of Green's function poles. Besides the self-energy can be expressed as 
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where Qi is the ith pole of the self-energy, is the weight of the pole, and Ns is the 
number of self-energy poles. 

Then the Dyson's equation (28) becomes 

y. 1 

S(oo) = + /i - eimp - 



First consider the S(oo), when a; — )■ oo 



r[iuj - Pi 



(A.3) 



S(cj — )■ oo) = iuj + ^ — e 



imp 



No y. 



iu + — ei 



imp 



Ng 



lUJ 



St 

i=i ^ 



Vi 



iu) 



Ng 



Ng 



i=l i=l 
Ng 



-1 



iu + fi- eimp - iu 

Ng 
i=l 



(A.4) 



For the the self energy poles in (A.3), it corresponding to the zero in the square 
bracket of the right side of (A.3), that is for Q G {Qi}-, we have 
Ng y 

The left hand side is monotonically decreasing between two adjacent Green's function 
poles, so we could use bisection method to find the self energy poles. 
We rewrite (A.3) as 



Ns 

E 



r[iu ~Qi 
For the left hand side 



iu + ji — e — S(oo) 



tr{iu-Pi 



(A.6) 



tl - Qj 



{iu - Qi 



lim 



' Na 



a y. 

lu + IJ,- eimp- T.{oo) - \J2 ^ 



'Ng y 

\im -{iu -Qi) ^ 



{iu - Qi) 



^^lU-Pj 

here we define 5 = iu — Qi and h{S) = Y.f=i q.^s-p ' 



(A,7) 



expand h{6) as 



(A.8) 
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so 



Ng 

E 



(Q, - P,) 

here we have got both self energy poles and its weight. 



(A.9) 
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